install.packages("mvtnorm")
library(mvtnorm)

install.packages("VGAM")
library(VGAM)

install.packages("trust")
library(trust)

library(MASS)

install.packages("ordinal")
library(ordinal)

library(Matrix)

library(mgcv)

install.packages("polycor")
library(polycor)


rm(list=ls(all=TRUE))

##############################

###########
# Dataset #
###########

setwd("C:\\Users\\PC\\Desktop\\PC PhD UCL, 31.08\\Thesis\\3 - Semi-parametric Polychotomous Ordinal Regression\\Application Drinking Behaviour")
#load("BCS70.drk1.RData")
load("BCS70.drk2.RData")

#DS <- BCS70.drk1
DS <- BCS70.drk2

setwd("C:\\Users\\PC\\Desktop\\PC PhD UCL, 31.08\\Thesis\\SemiParCLM\\10.17\\SemiParCLM - 10.17\\BIVOrdProbit")

source("1 - st.val.R")
source("2 - LkhoodOnly.R")
source("3 - Gradient.R")
source("4 - Hessian.R")
source("5 - LkhoodFull.R")
source("6 - BivCLM.R")


setwd("C:\\Users\\PC\\Desktop\\PC PhD UCL, 31.08\\Thesis\\SemiParCLM\\10.17\\SemiParCLM - 10.17")

source("01 - frm.sem.R")
source("02 - st.val.sem.ind.R")
source("02 - st.val.sem.R")
source("03 - lk.sem.np.R")
source("04 - pen.sem.R")
source("05 - gr.sem.D.R")
source("06 - gr.sem.u.R")
source("07 - hs.sem.W.R")
source("08 - hs.sem.K.R")
source("09 - lk.sem.p.R")
source("10 - lk.sem.nofit.R")
source("11 - lk.sem.fit.R")
source("12 - wlm.gen.R")
source("12 - wlm.R")
source("13 - no.conv.R")
source("14 - SemiParCLM.fit.R")
source("15 - pp.R")
source("16 - post.est.R")
source("17 - SemiParCLM.R")
source("18 - plot1.SemiParCLM.R")
source("19 - plot2.SemiParCLM.R")
source("20 - inference.sem.R")


##############
# Covariates #
##############

edu          <- DS$edu
drk          <- DS$drk
drk5         <- DS$drk5

mum.not.pres <- DS$mum.not.pres
dad.not.pres <- DS$dad.not.pres
mum.edu      <- DS$mum.edu
dad.edu      <- DS$dad.edu
s.class      <- DS$s.class
eth.child    <- DS$eth.child
mum.int.edu  <- DS$mum.int.edu
dad.int.edu  <- DS$dad.int.edu
sex.b        <- DS$sex.b
region       <- DS$region
home         <- DS$home

mum.wrk.hr <- DS$mum.wrk.hr ; BAS.tot <- DS$BAS.tot


############
# Analysis #
############

### Drinking behaviour ###

formula.1 <- edu  ~ mum.not.pres+dad.not.pres+s(mum.wrk.hr)+mum.edu+dad.edu+
                    s.class+s(BAS.tot)+eth.child+mum.int.edu+dad.int.edu+sex.b+home
                   
formula.2 <- drk5 ~ mum.not.pres+dad.not.pres+s(mum.wrk.hr)+mum.edu+dad.edu+
                    s.class+eth.child+mum.int.edu+dad.int.edu+sex.b+region+home
   
formula <- list(formula.1,formula.2)


ptm <- proc.time()

res.sem <- SemiParCLM(formula,iter.max=100,rint=1,rmax=5,poly=F,gcv=F,gamma=1.2,ind=F,wlm.g=F,semipar=T)

proc.time() - ptm

par     <- res.sem$par     ; value <- res.sem$value ; gradient  <- res.sem$gradient
hessian <- res.sem$hessian ; sp    <- res.sem$sp    ; converged <- res.sem$converged 
PP      <- res.sem$PP      ; se    <- res.sem$se

par(mfrow=c(1,3), mar=c(4,4.7,3.1,2), cex.lab = 1.6, cex.axis=1.3)
plot1.SemiParCLM(res.sem,shade=T,pages=0,scale=0,seWithMean=T)#,ylim=c(-1.5,2.2))
plot2.SemiParCLM(res.sem,shade=T,pages=0,scale=0,seWithMean=T)#,ylim=c(-1.5,2.2))

dev.copy(postscript, "smoothdrink2new.eps", height=3.7)
dev.off() 


#########################
# Analysis 2: Shrinkage #
#########################

### Smooth selection: shrinkage ###

formula.1 <- edu  ~ mum.not.pres+dad.not.pres+s(mum.wrk.hr,bs="ts")+mum.edu+dad.edu+
                    s.class+s(BAS.tot,bs="ts")+eth.child+mum.int.edu+dad.int.edu+sex.b+home                   
formula.2 <- drk5 ~ mum.not.pres+dad.not.pres+s(mum.wrk.hr,bs="ts")+mum.edu+dad.edu+
                    s.class+eth.child+mum.int.edu+dad.int.edu+sex.b+region+home
formula   <- list(formula.1,formula.2)

res.sem <- SemiParCLM(formula,iter.max=100,rint=1,rmax=5,poly=F,gcv=F,gamma=1.2,ind=F,wlm.g=F,semipar=T)

par(mfrow=c(1,3), mar=c(4,4.7,3.1,2), cex.lab = 1.6, cex.axis=1.3)
plot1.SemiParCLM(res.sem,shade=T,pages=0,scale=0,seWithMean=TRUE)
plot2.SemiParCLM(res.sem,shade=T,pages=0,scale=0,seWithMean=TRUE)

dev.copy(postscript, "Shrinkage.eps", height=3.7)
dev.off() 


### Drinking behaviour: one less smooth ###

formula.1 <- edu  ~ mum.not.pres+dad.not.pres+s(mum.wrk.hr)+mum.edu+dad.edu+
                    s.class+s(BAS.tot)+eth.child+mum.int.edu+dad.int.edu+sex.b+home
                   
formula.2 <- drk5 ~ mum.not.pres+dad.not.pres+mum.edu+dad.edu+
                    s.class+eth.child+mum.int.edu+dad.int.edu+sex.b+region+home
   
formula <- list(formula.1,formula.2)

#iter.max=100;rint=1;rmax=5;poly=F;gcv=F;gamma=1.2;ind=F;wlm.g=F;semipar=T
res.sem <- SemiParCLM(formula,iter.max=100,rint=1,rmax=5,poly=F,gcv=F,gamma=1.2,ind=F,wlm.g=F,semipar=T)

par     <- res.sem$par     ; value <- res.sem$value ; gradient  <- res.sem$gradient
hessian <- res.sem$hessian ; sp    <- res.sem$sp    ; converged <- res.sem$converged 
PP      <- res.sem$PP      ; se    <- res.sem$se

J.inv   <- res.sem$J.inv   ; V.inf <- res.sem$V.inf
K1      <- max(DS$edu)     ; K2 <- max(DS$drk5)

CI <- inference.sem(alpha=0.05,J.inv=J.inv,par=par,se=se,den=V.inf$den,K1=K1,K2=K2,lin.pred1=V.inf$lin.pred1,lin.pred2=V.inf$lin.pred2,rho.bar=V.inf$rho.bar)
CI.PP <- CI$CI.PP

plot1.SemiParCLM(res.sem,shade=T,pages=0,scale=0,seWithMean=T)#,ylim=c(-1.5,2.2))



